'simple fish frequency pdf deconvolution c.s. clay 5/4/93
'reads a file of fish echo data that is binned in logrithmic bins.
'echo frequencies are number per ping.

DIM WT(90),wF(90), bm(90),f(90),v(90),dv(90)
DIM ue(90),uF(90)
PI = 4 *  ATN (1):LD = 1 /  LOG (10): ka2 =10
da = LOG(10)/10
'default
r1 = 95 : r2 = 105 

100 PRINT" 't'    for transducer pdf"
PRINT" 'st'    for stored transducer pdf"
PRINT "'r'     for range gate and number of pings"
PRINT" 'f'    for fish frequency f(m) file and  compute u(m)"
PRINT" 'c'    to compute"
PRINT" 'mf'    to make file"
PRINT" 'q'    to quit"
INPUT q$

110 IF q$ = "t" GOTO 200
IF q$ = "f" GOTO 300
IF q$ = "r" GOTO 400
IF q$ = "c" GOTO 500
IF q$ = "mf" GOTO 600
IF q$ = "st" GOTO 800
IF q$ = "q" GOTO 1000
GOTO 100

200  PRINT " READ transducer FILE NAMED";:INPUT name1$
     
     OPEN name1$ FOR INPUT AS #1
     
INPUT #1, MD:
INPUT #1, ka2

FOR j = 0 TO MD
    INPUT #1, bm(j),W
    WT(j) = W/ka2^2
    PRINT j, WT(j)
NEXT j
CLOSE #1
INPUT q$
GOTO 100

300 PRINT " READ fish FILE NAMED";:INPUT name2$
'f(j) = number of echoes per ping
     
     OPEN name2$ FOR INPUT AS #2
     
INPUT #2, BJ:
INPUT #2, ka2
INPUT #2, r1,r2
sum = 0
FOR j = 0 TO BJ
    INPUT #2, v(j),f(j)
    sumf = sumf + f(j)    
    PRINT j,f(j)
NEXT j
PRINT"sumf =";sumf
CLOSE #2

vg = 2*PI*(r2^3-r1^3)/3

'compute the voltage bins for the echo frequency data

FOR m = 1 TO BJ+1
    dv(m) = EXP(-(m-.5)*da)-EXP(-(m+.5)*da)
NEXT m
ue(0)=0

'compute the distribution function for the echoes

FOR j = 1 TO BJ
     ue(j) = f(j)/(vg*dv(j))
     PRINT j, f(j), ue(j)
NEXT j

INPUT q$
GOTO 100

400 'new data
PRINT "old r1 and r2 =";r1;",";r2;" new = ";:INPUT r1,r2
GOTO 100

500 ' DECONVOLUTION
 
'  DECONVOLVE WE(E) TO GET WF(E)
'  USE RECURSIVE POLYNOMIAL LONG DIVISION ALGORITHM.
    PRINT 
    PRINT "DECONVOLVE WE(E)/WT(B)"
    'Since fish density Nf is not known, deconvolve uF(e) = uE(e)/wT(b)
    'Nf is computed after deconvolution is finished. 
    
    uF(0) = ue(0) / WT(0)
   FOR j = 1 TO BJ-1
       S = 0
          FOR m = 1 TO MD
              SC = 0
            IF j >  = m THEN SC = WT(m) * uF(j - m)
              S = S + SC
          NEXT m
          uF(j) =  - S / WT(0)
        IF j < MD THEN uF(j) = uF(j) + ue(j) / WT(0)
   NEXT j
   
'divide wF(j) by da = log(10)/nd and print

PRINT "j","uF(j)"
FOR j= 0 TO BJ-1
     uF(j) = uF(j)/da
PRINT j, uF(j)
NEXT j

PRINT"For density computation sum from j1 to j2 =";:INPUT j1,j2
 sumuf = 0
 FOR j = j1 TO j2
     sumuf = sumuf + uF(j)*dv(j)
 NEXT j
 
nf = sumuf :'nf = fish density

FOR j= 0 TO BJ-1
    wF(j) = uF(j)/nf
 PRINT wF(j)
NEXT j

PRINT "using wF(j), fish density =";nf

INPUT q$
PRINT
GOTO 100


600 'make file
PRINT " name THE FILE";:INPUT name3$

OPEN name3$ FOR OUTPUT AS #3
WRITE #3,name1$
WRITE #3,name2$
WRITE #3,BJ
WRITE #3,"m,  b(m),   wF(m)"
  FOR m = 0 TO BJ
    WRITE #3,bm(m), wF(m)
  NEXT m
WRITE #3,"r1, r2, mp =",r1,r2
WRITE #3,"fish density",nf
WRITE #3,"m,  ue(m),   dv(m)"
FOR m = 0 TO BJ
    WRITE #3,m, ue(m), dv(m)
NEXT m
CLOSE #3
 GOTO 100
 
 800 'transducer data file for wt(b) 2 significant figures (ka = 10)
name1$ = "stored wt(m)"
 PRINT" old ka =";ka2;"  new ka =";:INPUT ka2  
 MD=11:bm(0)=1:bm(1)=.79:bm(2)=.63:bm(3)=.5:bm(4)=.4:bm(5)=.32:bm(6)=.25
     bm(7)=.2:bm(8)=.16:bm(9)=.13:bm(10)=.1:bm(11)=.079
     WT(0)=2!:WT(1)=2.4:WT(2)=3!:WT(3)=3.6:WT(4)=4.3:WT(5)=5.2:WT(6)=6.3
     WT(7)=7.6:WT(8)=9.1:WT(9)=11:WT(10)=13:WT(11)=15
     
     FOR j = 0 TO MD
        WT(j) = WT(j)/ka2^2
     NEXT j
     
 GOTO 100
 
1000 END
